# Set the working directory

library(magrittr)
library(tidyverse)
library(dplyr)
library(scales)
library(ggplot2)
library(LearnBayes)
library(dampack)
source("~/Dropbox/Accueil Pasteur/Hépatite B/treat B/outputs/R code/PSA/DALY_function_PSA.R")
source("~/Dropbox/Accueil Pasteur/Hépatite B/treat B/outputs/R code/PSA/cost_functions_PSA.R")
source("~/Dropbox/Accueil Pasteur/Hépatite B/treat B/outputs/R code/PSA/filterint_incidence CC_HCC_PSA.R")
wd_root <- "~/Dropbox/Accueil Pasteur/Hépatite B/treat B/outputs/ultimate analysis/PSA/runs"
setwd(wd_root)


list_folder <- list.dirs(wd_root, recursive = FALSE)

#read CSV
Main_trace <- list()
Name_trace <- c()
for (i in c(1:length(list_folder))){
  wd <-  paste(list_folder[i],"/Main_trace",sep = "")
  print(wd)
  Main_trace[[i]] <- get_markov(wd)
  Name_trace[i] <- substr(list_folder[i],100,(nchar(list_folder[i])))
}
names(Main_trace) <- Name_trace


#get Years of life saved
list_YLS <- Main_trace
names_iteration <- names(Main_trace[[1]])
names_HS <- names(list_YLS$TREATB$`1`)

for (i in names(list_YLS)){
  for (j in names_iteration){
    a <- as.data.frame(Main_trace[[i]][j]) %>%
      setNames(names_HS)
    b <- as.data.frame(Main_trace[["WHO"]][j])%>%
      setNames(names_HS)
    list_YLS[[i]][j] <- get_YLS(b,a)
  }
}

list_YLS_clean <- lapply(list_YLS,function(x){
  sapply(
    indices_to_keep,function(i){
      x[[i]] 
    })
})
YLS_mean <- list()
YLS_LB <- list()
YLS_UB <- list()
for (i in names(list_YLS_clean)){
  YLS_mean[i] <- unlist(list_YLS_clean[[i]]) %>%
      mean()
  YLS_LB[i] <- unlist(list_YLS_clean[[i]]) %>%
    quantile(probs = c(0.025))
  YLS_UB[i] <- unlist(list_YLS_clean[[i]]) %>%
    quantile(probs = c(0.975))
}

YLS_df <- as.data.frame(do.call(cbind,list_YLS_clean), row.names = rownames(list_YLS_clean))
colnames(YLS_df) <- c("EASL","NatHist","Treat all","TREAT-B","WHO")
YLS_long <- YLS_df %>% 
  pivot_longer(cols = everything(), names_to = "Group", values_to = "Value") %>%
  filter(Group != "WHO") %>%
  filter(Group != "NatHist") %>%
  mutate(Value = as.numeric(Value)) %>%
  filter(!is.na(Value))

ggplot(YLS_long, aes(x = Value, fill = Group)) +
  geom_density(alpha = 0.5) +
  labs(title = "Density Plot of YLS", x = "Value", y = "Density") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black")+
  theme(text = element_text(size = 20)) 

YLS_it <- cbind(YLS_df,"iteration" = c(1:length(indices_to_keep)))
YLS_it_long <- YLS_it %>% 
  pivot_longer(cols = c("EASL", "TREAT-B","Treat all","WHO","NatHist"), names_to = "Group", values_to = "Value") %>%
  filter(Group != "WHO") %>%
  filter(Group != "NatHist") %>%
  mutate(Value = as.numeric(Value)) %>%
  filter(!is.na(Value))

ggplot(YLS_it_long, aes(x = iteration, y= Value, color= Group )) +
  geom_line() +
  labs(title = "YLS by group", x = "iteration", y = "YLS") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black")+
  theme(text = element_text(size = 20)) 

#get Years of life saved, discounted
list_YLS_disc <- Main_trace

for (i in names(list_YLS_disc)){
  for (j in names_iteration){
    list_YLS_disc[[i]][j] <- get_YLS_disc(Markov_BC = list_YLS_disc[["WHO"]][[j]],Markov_strat = list_YLS_disc[[i]][[j]])
  }
}

# remove iteration with out of range calibration targets

list_YLS_disc_clean <- lapply(list_YLS_disc,function(x){
  sapply(
    indices_to_keep,function(i){
      x[[i]] 
    })
})
colnames(YLS_df) <- c("EASL","NatHist","Treat all","TREAT-B","WHO")  

YLS_disc_mean <- list()
YLS_disc_LB <- list()
YLS_disc_UB <- list()
for (i in names(list_YLS_disc_clean)){
  YLS_disc_mean[i] <- unlist(list_YLS_disc_clean[[i]]) %>%
    mean()
  YLS_disc_LB[i] <- unlist(list_YLS_disc_clean[[i]]) %>%
    quantile(probs = c(0.025))
  YLS_disc_UB[i] <- unlist(list_YLS_disc_clean[[i]]) %>%
    quantile(probs = c(0.975))
}

YLS_disc_df <- as.data.frame(do.call(cbind,list_YLS_disc_clean), row.names = rownames(list_YLS_disc_clean))
colnames(YLS_disc_df) <- c("EASL","NatHist","Treat all","TREAT-B","WHO")  
YLS_disc_long <- YLS_disc_df %>% 
  pivot_longer(cols = everything(), names_to = "Group", values_to = "Value") %>%
  filter(Group != "WHO") %>%
  filter(Group != "NatHist") %>%
  mutate(Value = as.numeric(Value)) %>%
  filter(!is.na(Value))

ggplot(YLS_disc_long, aes(x = Value, fill = Group)) +
  geom_density(alpha = 0.5) +
  labs(title = "Density Plot of Years of Life Saved, discounted", x = "Value", y = "Density") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black")+
  scale_y_continuous(labels = function(x) format(x, scientific = TRUE))+
  theme(text = element_text(size = 12)) 


YLS_disc_it <- cbind(YLS_disc_df,"iteration" = c(1:length(indices_to_keep)))
YLS_disc_it_long <- YLS_disc_it %>% 
  pivot_longer(cols = c("EASL", "TREAT-B","Treat all","WHO","NatHist"), names_to = "Group", values_to = "Value") %>%
  filter(Group != "WHO") %>%
  filter(Group != "NatHist") %>%
  mutate(Value = as.numeric(Value)) %>%
  filter(!is.na(Value))


ggplot(YLS_disc_it_long, aes(x = iteration, y= Value, color= Group )) +
  geom_line() +
  labs(title = "YLS disc by group", x = "iteration", y = "YLS") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black")+
  theme(text = element_text(size = 20)) 


# Create disability Weight index

#find the beta distribution of disability weight using the root-finding algorithm

#DCC
DCC_mean <- 0.178
DCC_lower_ci <- 0.123
DCC_upper_ci <- 0.259

quantile1=list(p=.025, x=DCC_lower_ci)    
quantile2=list(p=.975, x=DCC_upper_ci)  
DCC_bdistrib <- beta.select(quantile1, quantile2)
DCC_PSA_value <- rbeta(500, DCC_bdistrib[1], DCC_bdistrib[2])

#HCC
HCC_mean <- 0.288
HCC_lower_ci <- 0.193
HCC_upper_ci <- 0.399

quantile1=list(p=.025, x=HCC_lower_ci)  
quantile2=list(p=.975, x=HCC_upper_ci)  
HCC_bdistrib <- beta.select(quantile1, quantile2)
HCC_PSA_value <- rbeta(500, HCC_bdistrib[1], HCC_bdistrib[2])
##########################@
col_names <- names_HS
Dw_it_list <- list(rep(0,500))
DW_list_trt_gbd <- rep(Dw_it_list,37)
names(DW_list_trt_gbd) <- col_names[1:(length(col_names)-1)]

for (i in c(1:length(DW_list_trt_gbd))) {
  print(DW_list_trt_gbd[[i]])
  DW_list_trt_gbd[[i]] <- rep(0,500)
  if (grepl("\\.trt",names(DW_list_trt_gbd[i])) == TRUE) DW_list_trt_gbd[[i]] <- rep(0,500)
  if (grepl("CC.trt",names(DW_list_trt_gbd[i])) == TRUE) DW_list_trt_gbd[[i]] <- rep(0,500)
  if (grepl("DCC",names(DW_list_trt_gbd[i])) == TRUE) DW_list_trt_gbd[[i]] <- DCC_PSA_value
  if (grepl("HCC",names(DW_list_trt_gbd[i])) == TRUE) DW_list_trt_gbd[[i]] <- HCC_PSA_value
}

####COST####

# Calculate cost

# PSA using gamma distribution with hypothesis of sd= 20% mean

#admission cost
mean_hospit_admission_cost <- 74.7
sd_hospit_admission_cost <- 0.3*(mean_hospit_admission_cost)
theta_hospit_admission_cost <- (sd_hospit_admission_cost^2)/mean_hospit_admission_cost
kappa_H_hospit_admission_cost <- mean_hospit_admission_cost/theta_hospit_admission_cost

hospit_admission_cost <- rgamma(500,shape=kappa_H_hospit_admission_cost,scale=theta_hospit_admission_cost)

#nb of admission per year
nb_hospi_CC <- runif(500,0,4)
nb_hospi_DCC <- runif(500,0,6)
nb_hospi_HCC <- runif(500,0,6)

CC_cost <- hospit_admission_cost*nb_hospi_CC
DCC_cost <-  hospit_admission_cost*nb_hospi_DCC
HCC_cost <-  hospit_admission_cost*nb_hospi_HCC
#Cost on treatment

mean_treat_cost <- 55.1
sd_treat_cost <- 0.2*(mean_treat_cost)
theta_treat_cost <- (sd_treat_cost^2)/mean_treat_cost
kappa_H_treat_cost <- mean_treat_cost/theta_treat_cost
treat_cost <- rgamma(500,shape=kappa_H_treat_cost,scale=theta_treat_cost)

#TREATB HBe Positif, off treatment

mean_HBePos_TREATB_cost <- 21.7
sd_HBePos_TREATB_cost <- 0.2*(mean_HBePos_TREATB_cost)
theta_HBePos_TREATB_cost <- (sd_HBePos_TREATB_cost^2)/mean_HBePos_TREATB_cost
kappa_H_HBePos_TREATB_cost <- mean_HBePos_TREATB_cost/theta_HBePos_TREATB_cost
HBePos_TREATB_cost <- rgamma(500,shape=kappa_H_HBePos_TREATB_cost,scale=theta_HBePos_TREATB_cost)

#TREATB HBe neg, off treatment

mean_HBeNeg_TREATB_cost <- 14.2
sd_HBeNeg_TREATB_cost <- 0.2*(mean_HBeNeg_TREATB_cost)
theta_HBeNeg_TREATB_cost <- (sd_HBeNeg_TREATB_cost^2)/mean_HBeNeg_TREATB_cost
kappa_H_HBeNeg_TREATB_cost <- mean_HBeNeg_TREATB_cost/theta_HBeNeg_TREATB_cost
HBeNeg_TREATB_cost <- rgamma(500,shape=kappa_H_HBeNeg_TREATB_cost,scale=theta_HBeNeg_TREATB_cost)

#WHO HBe Positif, off treatment

mean_HBePos_WHO_cost <- 53.0
sd_HBePos_WHO_cost <- 0.2*(mean_HBePos_WHO_cost)
theta_HBePos_WHO_cost <- (sd_HBePos_WHO_cost^2)/mean_HBePos_WHO_cost
kappa_H_HBePos_WHO_cost <- mean_HBePos_WHO_cost/theta_HBePos_WHO_cost
HBePos_WHO_cost <- rgamma(500,shape=kappa_H_HBePos_WHO_cost,scale=theta_HBePos_WHO_cost)

#WHO HBe neg, off treatment

mean_HBeNeg_WHO_cost <- 45.5
sd_HBeNeg_WHO_cost <- 0.2*(mean_HBeNeg_WHO_cost)
theta_HBeNeg_WHO_cost <- (sd_HBeNeg_WHO_cost^2)/mean_HBeNeg_WHO_cost
kappa_H_HBeNeg_WHO_cost <- mean_HBeNeg_WHO_cost/theta_HBeNeg_WHO_cost
HBeNeg_WHO_cost <- rgamma(500,shape=kappa_H_HBeNeg_WHO_cost,scale=theta_HBeNeg_WHO_cost)

#EASL HBe Positif, off treatment

mean_HBePos_EASL_cost <- 59.1
sd_HBePos_EASL_cost <- 0.2*(mean_HBePos_EASL_cost)
theta_HBePos_EASL_cost <- (sd_HBePos_EASL_cost^2)/mean_HBePos_EASL_cost
kappa_H_HBePos_EASL_cost <- mean_HBePos_EASL_cost/theta_HBePos_EASL_cost
HBePos_EASL_cost <- rgamma(500,shape=kappa_H_HBePos_EASL_cost,scale=theta_HBePos_EASL_cost)

#EASL HBe neg, off treatment

mean_HBeNeg_EASL_cost <- 51.6
sd_HBeNeg_EASL_cost <- 0.2*(mean_HBeNeg_EASL_cost)
theta_HBeNeg_EASL_cost <- (sd_HBeNeg_EASL_cost^2)/mean_HBeNeg_EASL_cost
kappa_H_HBeNeg_EASL_cost <- mean_HBeNeg_EASL_cost/theta_HBeNeg_EASL_cost
HBeNeg_EASL_cost <- rgamma(500,shape=kappa_H_HBeNeg_EASL_cost,scale=theta_HBeNeg_EASL_cost)


#created nested list of costs
it_list <- list(rep(0,500))
l_cost <- rep(it_list,7)
cost_names <- c("out_ntrt_EP","out_ntrt_EN","out_trt_EP", "out_trt_EN","Cirr","DCC","HCC")
names(l_cost) <- cost_names

l_v_cost_NH <- l_cost %>%
  modify_at("Cirr", ~CC_cost) %>%
  modify_at("DCC", ~DCC_cost) %>%
  modify_at("HCC", ~HCC_cost)
  

l_v_cost_TREATB <- l_cost %>%
  modify_at("out_ntrt_EP", ~HBePos_TREATB_cost) %>%
  modify_at("out_ntrt_EN", ~HBeNeg_TREATB_cost) %>%
  modify_at("out_trt_EP", ~treat_cost) %>%
  modify_at("out_trt_EN", ~treat_cost) %>%
  modify_at("Cirr", ~CC_cost) %>%
  modify_at("DCC", ~DCC_cost) %>%
  modify_at("HCC", ~HCC_cost)

l_v_cost_WHO <- l_cost %>%
  modify_at("out_ntrt_EP", ~HBePos_WHO_cost) %>%
  modify_at("out_ntrt_EN", ~HBeNeg_WHO_cost) %>%
  modify_at("out_trt_EP", ~treat_cost) %>%
  modify_at("out_trt_EN", ~treat_cost) %>%
  modify_at("Cirr", ~CC_cost) %>%
  modify_at("DCC", ~DCC_cost) %>%
  modify_at("HCC", ~HCC_cost)

l_v_cost_EASL <- l_cost %>%
  modify_at("out_ntrt_EP", ~HBePos_EASL_cost) %>%
  modify_at("out_ntrt_EN", ~HBeNeg_EASL_cost) %>%
  modify_at("out_trt_EP", ~treat_cost) %>%
  modify_at("out_trt_EN", ~treat_cost) %>%
  modify_at("Cirr", ~CC_cost) %>%
  modify_at("DCC", ~DCC_cost) %>%
  modify_at("HCC", ~HCC_cost)

l_v_cost_TREATALL <- l_cost %>%
  modify_at("out_trt_EP", ~treat_cost) %>%
  modify_at("out_trt_EN", ~treat_cost) %>%
  modify_at("Cirr", ~CC_cost) %>%
  modify_at("DCC", ~DCC_cost) %>%
  modify_at("HCC", ~HCC_cost)

cost_it_list <- list(rep(0,500))
cost_list_main_trace <- rep(cost_it_list,37)
names(cost_list_main_trace) <- names_HS[1:(length(col_names)-1)]

l_cost_NatHist <- get_l_variable_cost(cost_list_main_trace, l_v_cost_NH)
l_cost_TREATB <- get_l_variable_cost(cost_list_main_trace,l_v_cost_TREATB)
l_cost_WHO <- get_l_variable_cost(cost_list_main_trace,l_v_cost_WHO)
l_cost_EASL <- get_l_variable_cost(cost_list_main_trace,l_v_cost_EASL)
l_cost_TREATALL <- get_l_variable_cost(cost_list_main_trace,l_v_cost_TREATALL)

#variable cost discount


get_var_cost_disc_strat <- function(strat){
  Trace <- Main_trace[[strat]]
  l_cost_strat <- get(paste0("l_cost_",strat))
  var_cost_disc <- list()
  for (i in names(Trace)){
    df_cost <- map(l_cost_strat,~.x [as.numeric(i)])
    var_cost_disc[i] <- get_var_cost_disc(as.data.frame(Trace[[i]]),df_cost)
  }
  return(var_cost_disc[as.character(1:500)])
}

var_cost_disc_EASL <- get_var_cost_disc_strat("EASL")
var_cost_disc_NatHist <- get_var_cost_disc_strat("NatHist")
var_cost_disc_WHO <- get_var_cost_disc_strat("WHO")
var_cost_disc_TREATALL <- get_var_cost_disc_strat("TREATALL")
var_cost_disc_TREATB <- get_var_cost_disc_strat("TREATB")

#Old fonction
"Trace_NH <- Main_trace$NatHist
var_cost_disc_NH <- list()
for (i in names(Trace_NH)){
  print(i)
  df_cost_NH <- map(l_cost_NH,~.x [as.numeric(i)])
  print(df_cost_NH)
  var_cost_disc_NH[i] <- get_var_cost_disc(as.data.frame(Trace_NH[[i]]),df_cost_NH)
}
"

#Graph representation
variable_cost_df <- as.data.frame(cbind("EASL" =var_cost_disc_EASL, "WHO" = var_cost_disc_WHO,
                                        "TREATB" = var_cost_disc_TREATB,"TREATALL"=var_cost_disc_TREATALL))
l_variable_cost <-  list("EASL" =var_cost_disc_EASL, "WHO" = var_cost_disc_WHO,
                             "TREATB" = var_cost_disc_TREATB,"TREATALL"=var_cost_disc_TREATALL)


variable_cost_df_long <- variable_cost_df %>% 
  pivot_longer(cols = everything(), names_to = "Group", values_to = "Value") %>%
  mutate(Value = as.numeric(Value)) %>%
  filter(!is.na(Value))

ggplot(variable_cost_df_long, aes(x = Value, fill = Group)) +
  geom_density(alpha = 0.5) +
  labs(title = "Density Plot of variable cost", x = "Value", y = "Density") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black")

variable_cost_df_it <- cbind(variable_cost_df,"iteration" = c(1:500))
variable_cost_df_it_long <- variable_cost_df_it %>% 
  pivot_longer(cols = c("EASL", "TREATB","TREATALL","WHO"), names_to = "Group", values_to = "Value") %>%
  mutate(Value = as.numeric(Value)) %>%
  filter(!is.na(Value))

ggplot(variable_cost_df_it_long, aes(x = iteration, y= Value, color= Group )) +
  geom_point() +
  labs(title = "variable cost by group", x = "iteration", y = "YLS") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black")

# implementing total costs
WHO_startup <- 0
WHO_fixed_cost <- 58072
EASL_startup <- 529760
EASL_fixed_cost <- 59752
discount <- exp(-0.03*(c(0:70)))
WHO_fixed_cost_disc <- sum(WHO_fixed_cost*discount)
EASL_fixed_cost_disc <- sum(EASL_fixed_cost*discount)

EASL_total_cost <- map(var_cost_disc_EASL, ~ .x + EASL_fixed_cost_disc)
WHO_total_cost <- map(var_cost_disc_WHO, ~ .x + EASL_fixed_cost_disc)

total_cost_df <- as.data.frame(cbind("EASL" = EASL_total_cost, "WHO" = WHO_total_cost,
                                       "TREATB" = var_cost_disc_TREATB,"TREATALL"=var_cost_disc_TREATALL))

total_incremental_cost_df <- lapply(total_cost_df, function(col) unlist(col) - unlist(total_cost_df$WHO))

l_total_cost <- list("EASL" = EASL_total_cost, "WHO" = WHO_total_cost,
                     "TREATB" = var_cost_disc_TREATB,"TREATALL"=var_cost_disc_TREATALL)

list_total_cost_clean <- lapply(l_total_cost,function(x){
  sapply(
    indices_to_keep,function(i){  x[[i]] })
})

l_total_incremental_cost <- lapply(l_total_cost, function(x) unlist(x) - unlist(WHO_total_cost))

list_incremental_cost_clean <- lapply(l_total_incremental_cost,function(x){
  sapply(
    indices_to_keep,function(i){  x[[i]] })
})

TREATB_incremental_cost <- list_incremental_cost_clean$TREATB
length(keep(TREATB_incremental_cost,TREATB_incremental_cost<0))
length(keep(list_incremental_cost_clean$EASL,list_incremental_cost_clean$EASL<0))/500
length(keep(list_incremental_cost_clean$TREATALL,list_incremental_cost_clean$TREATALL<0))
22total_cost_df_long <- as.data.frame(list_total_cost_clean) %>% 
  pivot_longer(cols = everything(), names_to = "Group", values_to = "Value") %>%
  mutate(Value = as.numeric(Value)) %>%
  filter(!is.na(Value))

ggplot(total_cost_df_long, aes(x = Value, fill = Group)) +
  geom_density(alpha = 0.5) +
  labs(title = "Density Plot of total cost", x = "Value", y = "Density") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black")

total_cost_df_it <- cbind(as.data.frame(list_total_cost_clean),"iteration" = c(1:length(indices_to_keep)))
total_cost_df_it_long <- total_cost_df_it %>% 
  pivot_longer(cols = c("EASL", "TREATB","TREATALL","WHO"), names_to = "Group", values_to = "Value") %>%
  mutate(Value = as.numeric(Value)) %>%
  filter(!is.na(Value))

ggplot(total_cost_df_it_long, aes(x = iteration, y= Value, color= Group )) +
  geom_point() +
  labs(title = "variable cost by group", x = "iteration", y = "YLS") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black")

#plot incremental cost
names(l_total_incremental_cost) <-  c("EASL", "WHO", "TREAT-B","Treat all")
total_incr_cost_df_long <- as_tibble(l_total_incremental_cost, col.names = names(l_total_incremental_cost)) %>%
  pivot_longer(cols = everything(), names_to = "Group", values_to = "Value") %>%
  mutate(Value = as.numeric(Value)) %>%
  filter(!is.na(Value)) %>%
  filter(Group != "WHO")

ggplot(total_incr_cost_df_long, aes(x = Value, color = Group)) +
  stat_ecdf() +
  #scale_x_continuous(limits = c(-100,100)) +
  labs(title = "Density Plot of incremental cost", x = "Value", y = "Density")+
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
  theme(text = element_text(size = 12))   

# Calculation of YLL

# YLL ---------------------------------------------------------------------


Main_trace_YLL <- lapply(Main_trace,function(x){
  sapply(
    names_iteration,function(i){
      return(get_YLL(x[[i]],LE = 61,age_init = 20))
    })
})

#calculate YLL, discounted


Main_trace_YLL_disc <- lapply(Main_trace,function(x){
  sapply(
    names_iteration,function(i){
      return(get_YLL_disc(x[[i]],LE = 61,age_init = 20))
    })
})

#Old function
#list_YLL_disc <- Main_trace
#for (i in c(1:length(Main_trace))){
#  for (j in c(1:length(Main_trace[[i]]))){
#    print(names(Main_trace[i]))
#    #print(names(Main_trace[[i]][j]))
#    list_YLL_disc[[i]][j] <- get_YLL(as.data.frame(Main_trace[[i]][j]) %>% setNames(names_HS),LE = 61,age_init = 20)
#  }
#}

#Calculation of YLD
Main_trace_YLD <- lapply(Main_trace,function(x){
  sapply(
    names_iteration,function(i){
      DW_list <- sapply(DW_list_trt_gbd,function(j){
        j[as.numeric(i)]
      })
      return(get_YLD(x[[i]],DW_list))
    })
})

#Calculation of YLD, discounted

Main_trace_YLD_disc <- lapply(Main_trace,function(x){
  sapply(
    names_iteration,function(i){
      DW_list <- sapply(DW_list_trt_gbd,function(j){
        j[as.numeric(i)]
      })
      return(get_YLD_disc(x[[i]],DW_list))
    })
})

# calculation of DALY

Main_trace_DALY <- map2(Main_trace_YLD, Main_trace_YLL, ~ .x + .y)
Main_trace_DALY_disc <- map2(Main_trace_YLD_disc, Main_trace_YLL_disc, ~ .x + .y)
Main_trace_DALY_disc_clean <- lapply(Main_trace_DALY_disc,function(x){
  sapply(
    indices_to_keep,function(i){  x[[i]] })
})
Main_trace_DALY_averted_disc <- map(Main_trace_DALY_disc, ~ Main_trace_DALY_disc$WHO -.x )
Main_trace_DALY_averted_disc_clean <- lapply(Main_trace_DALY_averted_disc,function(x){
  sapply(
    indices_to_keep,function(i){  x[[i]] })
})

mean(ICER_TREATALL)*1.12
quantile(ICER_TREATALL, 0.975)*1.12
quantile(ICER_TREATALL, 0.025)*1.12
#calculation of ICER


Main_trace_ICER <- Main_trace_DALY_averted_disc_clean  
Main_trace_ICER$EASL <- map2(Main_trace_DALY_averted_disc_clean$EASL, list_incremental_cost_clean$EASL , ~ .y / .x) 
Main_trace_ICER$TREATALL <- map2(Main_trace_DALY_averted_disc_clean$TREATALL, list_incremental_cost_clean$TREATALL, ~ .y / .x) 
Main_trace_ICER$TREATB <- map2(Main_trace_DALY_averted_disc_clean$TREATB,list_incremental_cost_clean$TREATB , ~ .y / .x) 
Main_trace_ICER$NatHist <- NULL
Main_trace_ICER$WHO <- NULL

Main_trace_ICER_YLS <- list_YLS_disc_clean
Main_trace_ICER_YLS$EASL <- map2(list_YLS_disc_clean$EASL, list_incremental_cost_clean$EASL , ~ .y / .x) 
Main_trace_ICER_YLS$TREATALL <- map2(list_YLS_disc_clean$TREATALL, list_incremental_cost_clean$TREATALL, ~ .y / .x) 
Main_trace_ICER_YLS$TREATB <- map2(list_YLS_disc_clean$TREATB,list_incremental_cost_clean$TREATB , ~ .y / .x) 
Main_trace_ICER_YLS$NatHist <- NULL
Main_trace_ICER_YLS$WHO <- NULL


ICER_df_long <- as.data.frame(do.call(cbind,Main_trace_ICER)) %>% 
  pivot_longer(cols = everything(), names_to = "Group", values_to = "Value") %>%
  mutate(Value = as.numeric(Value)) %>%
  filter(!is.na(Value))

ggplot(ICER_df_long, aes(x = Value, fill = Group)) +
  geom_density(alpha = 0.5) +
  scale_x_continuous(limits = c(-2000,2000)) +
  labs(title = "Density Plot of ICER", x = "Value", y = "Density")


plot(density(Main_trace_DALY_averted_disc_clean$TREATALL))
DALY_df <- as.data.frame(Main_trace_DALY_disc)

DALY_av_df <- as.data.frame(Main_trace_DALY_averted_disc_clean)
colnames(DALY_av_df) <- c("EASL","NatHist","Treat all","TREAT-B","WHO")  
DALY_av_df_long <- DALY_av_df %>% 
  pivot_longer(cols = everything(), names_to = "Group", values_to = "Value") %>%
  mutate(Value = as.numeric(Value)) %>%
  filter(!is.na(Value)) %>%
  filter(Group != "WHO") %>%
  filter(Group != "NatHist")

ggplot(DALY_av_df_long, aes(x = Value, fill = Group)) +
  geom_density(alpha = 0.5) +
  labs(title = "Density Plot of DALY averted", x = "Value", y = "Density") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black")+
  theme(text = element_text(size = 12)) 

DALY_av_df_it <- cbind(as.data.frame(Main_trace_DALY_averted_disc_clean),"iteration" = c(1:length(indices_to_keep)))
DALY_av_df_it_long <- DALY_av_df_it %>% 
  pivot_longer(cols = c("EASL", "TREATB","TREATALL","WHO"), names_to = "Group", values_to = "Value") %>%
  mutate(Value = as.numeric(Value)) %>%
  filter(!is.na(Value)) %>%
  filter(Group != "WHO") %>%
  filter(Group != "NatHist")

ggplot(DALY_av_df_it_long, aes(x = iteration, y= Value, color= Group )) +
  geom_point() +
  labs(title = "DALY averted by group", x = "iteration", y = "YLS") +
  #scale_y_continuous(limits = c(3.25*10^5,4*10^5)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black")


TREATALL_DALY_av <- Main_trace_DALY_averted_disc_clean$TREATALL
indices_pos_DALY_TREATALL_av <- which(TREATALL_DALY_av >0)

TREATALL_YLS_av <- list_YLS_disc_clean$TREATALL
indices_pos_YLS_TREATALL_av <- which(TREATALL_YLS_av >0)

length(indices_pos_DALY_TREATALL_av)
indices_pos_COST_TREATALL_av <- which(total_incremental_cost_df$TREATALL >0)
length(indices_pos_COST_TREATALL_av)

TREATALL_YLS_av <- list_YLS_disc_clean$TREATALL
indices_pos_YLS_TREATALL_av <- which(TREATALL_YLS_av >0)
length(indices_pos_YLS_TREATALL_av)
indices_pos_COST_TREATALL_av <- which(total_incremental_cost_df$TREATALL >0)
length(indices_pos_COST_TREATALL_av)

#remove negative cost and negative ouput
TREATALL_DALY_averted_disc <- Main_trace_DALY_averted_disc_clean$TREATALL[indices_pos_DALY_TREATALL_av]
TREATALL_DALY_cost_disc <- l_total_cost$TREATALL[indices_pos_DALY_TREATALL_av]
TREATALL_DALY_increment_cost_disc <- list_incremental_cost_clean$TREATALL[indices_pos_DALY_TREATALL_av]
TREATALL_DALY_increment_cost_disc_df <- as.data.frame(TREATALL_DALY_increment_cost_disc)
ggplot(TREATALL_DALY_increment_cost_disc_df, aes(x = TREATALL_DALY_increment_cost_disc)) +
  stat_ecdf() +
  #scale_x_continuous(limits = c(-100,100)) +
  labs(title = "Density Plot of treat all incremental cost", x = "Value", y = "Density")

ICER_TREATALL_av <- Main_trace_ICER$TREATALL[indices_pos_DALY_TREATALL_av]
ICER_TREATALL_av_YLS <- Main_trace_ICER_YLS$TREATALL[indices_pos_YLS_TREATALL_av]

ICER_TREATALL <- as.data.frame(do.call(cbind,ICER_TREATALL_av))
ICER_TREATALL <- unlist(ICER_TREATALL)
mean(ICER_TREATALL)
quantile(ICER_TREATALL, 0.975)
quantile(ICER_TREATALL, 0.025)

mean(unlist(ICER_TREATALL_av_YLS))
quantile(unlist(ICER_TREATALL_av_YLS), 0.975)
quantile(unlist(ICER_TREATALL_av_YLS), 0.025)
median(ICER_TREATALL)
ICER_TREATALL_df <-as.data.frame(ICER_TREATALL, row.names = as.character(c(1:length(ICER_TREATALL_av))))

ggplot(ICER_TREATALL_df, aes(x = ICER_TREATALL)) +
  #stat_ecdf() +
  geom_density() +
  scale_x_continuous(limits = c(-10000,50000)) +
  labs(title = "Cumulative incidence of TREAT ALL ICER", x = "Value", y = "Density")


wtp <- seq(from = 0, to = 200000, by= 500)


psa_obj_DALY_averted <- make_psa_obj(cost= as.data.frame(cbind("EASL" = unlist(list_incremental_cost_clean$EASL),
                                                               "Treat all" = unlist(list_incremental_cost_clean$TREATALL),
                                                               "TREAT-B" = unlist(list_incremental_cost_clean$TREATB))),
                                     effectiveness = as.data.frame(cbind("EASL" = unlist(Main_trace_DALY_averted_disc_clean$EASL),
                                                                         "Treat all" = unlist(Main_trace_DALY_averted_disc_clean$TREATALL),
                                                                         "TREAT-B" = unlist(Main_trace_DALY_averted_disc_clean$TREATB))),
                                     strategies = c("EASL","Treat all","TREAT-B")
)
plot(psa_obj_DALY_averted,txtsize = 12)

ceac_obj_DALY_averted <- ceac(wtp = wtp,psa = psa_obj_DALY_averted)
plot(ceac_obj_DALY_averted, 
     frontier = F, 
     points = T,txtsize = 12)

psa_obj_YLS_averted <- make_psa_obj(cost= as.data.frame(cbind("EASL" = unlist(list_incremental_cost_clean$EASL),
                                                              "TREATALL" = unlist(list_incremental_cost_clean$TREATALL),
                                                              "TREATB" = unlist(list_incremental_cost_clean$TREATB))),
                                    effectiveness = as.data.frame(cbind("EASL" = unlist(list_YLS_disc_clean$EASL),
                                                                        "TREATALL" = unlist(list_YLS_disc_clean$TREATALL),
                                                                        "TREATB" = unlist(list_YLS_disc_clean$TREATB))),
                                    strategies = c("EASL","TREATALL","TREATB")
)
plot(psa_obj_YLS_averted,txtsize = 12) 
calculate_icers(cost = psa_obj_YLS_averted$cost,effect = psa_obj_YLS_averted$effectiveness,strategies = psa_obj_YLS_averted$strategies)
mean(unlist(Main_trace_DALY_averted_disc_clean$TREATALL))
quantile(unlist(Main_trace_DALY_averted_disc_clean$TREATALL), 0.025)
quantile(unlist(Main_trace_DALY_averted_disc_clean$TREATALL), 0.975)
ceac_obj_YLS <- ceac(wtp = wtp,psa = psa_obj_YLS_averted)
plot(ceac_obj_YLS, 
     frontier = F, 
     points = TRUE)


